Linear prediction is well known for processing of signals in a number of applications. Using this method, the p+1th sample of a signal is predicted by forming a linear combination of p previous samples of the signal. The linear combination is usually optimized by minimizing the square of the prediction error. One of the most widely used models for linear prediction is the autoregressive (AR) model.
Perhaps the most widely recognized application of linear prediction is in the area of audio or speech coding, where an incoming signal stream is digitized, then samples of the digitized signal are processed by a variable linear predictive filter to determine relationships in the signal in order to calculate, i.e., predict, the value of the next output sample. This process of identifying model parameters corresponding to the signal permits the data in the signal to be compressed for faster transmission and also provides means for removing a significant amount of noise from the original received signal. Upon receipt of the transmission, the sound is regenerated using a sound synthesizer, which uses an excitation signal and information included in the transmission about the filter coefficients to restore the audio signal less any suppressed noise. Linear prediction has also been used in image processing and modeling of complex systems, including stock market analysis.
Linear prediction can be analyzed in the time or frequency domain. An important application of linear prediction is as a method of spectrum estimation. Generally, spectrum estimation techniques are used to extract useful data from a signal which also contains noise and/or interference. Specifically, spectrum estimation is a problem that involves estimating the power spectrum from a finite number of noisy measurements. “Useful data” includes not only information that may be encoded in the signal, but also characteristics of the signal itself, such as direction and delay, permitting signal detection and tracking. Thus, “useful data” is any information that is determinable from the signal. Such techniques also may be used for interpolation or extrapolation of missing portions of a signal. One exemplary application is in the area of wireless communication systems, where co-channel interference can be severe. Techniques which have been employed for increasing the signal level of the primary communication path include directive antennas at base station sites and adaptive beamforming, in which a phased antenna array is used to account for various angles between the base station and a remote user. While beamforming to a fixed location over a line of sight may be performed with relative ease, the task of transmitting to a mobile user over a time-varying multipath is more difficult. An adaptive transmit beamforming approach determines the angle of departure (AOD) at which energy is to be transmitted for the base station antenna array to a given remote user. Each AOD corresponds to one of the signal paths of the multipath channel, and is determined by estimating each angle of arrival (AOA) at the base station of signal energy from the user using a spectrum estimation technique. A transmit beam form is adaptively formed to maximize the signal level along each desired AOD while minimizing the radiation projected at all other angles. Several well known types of high-resolution or “super-resolution” methods, e.g., MUSIC, ESPRIT, Maximum likelihood and WSF, may be used to estimate an AOA spectrum in the presence of several signal sources for purposes of determining AOD.
Spectrum estimation has also been used for estimating time of arrival (TOA) which has applications in geolocation, e.g., triangulation. Super-resolution techniques such as MUSIC and root MUSIC have been used to provide significant improvement in resolution for TOA estimates. (See, e.g., U.S. Pat. No. 5,890,068, “Wireless Location System”, the disclosure of which is incorporated herein by reference.) TOA alone or in combination with AOA can be used for bearing estimation for sonar and radar applications (which utilize beamforming methods) and in synthetic aperture radar (SAR) imaging applications. Spectral estimation is also widely used in analytical methods including chromatographic and spectrographic analysis of materials, where multiple channels of data are contained within a measured spectrum and regression models are used to extract relevant data from irrelevant data and noise.
Because there are so many spectrum estimators, it is useful to compare a potential estimator to a theoretical limit given by the Cramer-Rao bound (CRB). The CRB gives the lower bound on the variance of the estimated parameter for unbiased estimators. For a single frequency in white noise, the conventional discrete Fourier transform (DFT) can achieve the CRB. For multiple, closely spaced frequencies, the DFT can no longer resolve the frequencies and other techniques have been developed that are capable of superior resolution. One such technique is the aforementioned autoregressive (AR) technique.
AR spectrum estimation is based on modeling a process as the output of an all-pole filter ap whose input is drawn from unit variance white noise. This is equivalent to forming a linear combination of the past p data values to estimate a new data value for a pth order process where ap is the unknown weight vector which this algorithm determines, i.e.,
                              x          ⁡                      (            n            )                          =                  -                                    ∑                              k                =                1                            p                        ⁢                                                            a                  p                  *                                ⁡                                  (                  k                  )                                            ⁢                                                x                  ⁡                                      (                                          n                      -                      k                                        )                                                  .                                                                        (        1        )            Describing the problem in terms of spatial frequencies, which represent angles of arrival of plane waves impinging on a uniform line array, the estimation process can be written as−aPHXp=d,   (2)where Xp is a p×k matrix whose columns represent k different snapshots of the array and whose p row contain the measured values of the first p array elements from a length (p+1) array. (“Snapshot” refers to a sample of data taken in an instant in time, e.g., a window or frame of data, such that a plurality of snapshots can represent time series data.) The vector x is a 1×k row vector composed of the data at the (p+1)th array element for each of the snapshots. The problem can be thought of in terms of linear prediction where a linear combination of p elements is used to form an estimate for the (p+1)th array element. The full rank solution, which minimizes the squared error between the actual and predicted values, assuming k>p, isa=−(XpHpH)−1XpxH,   (3)which can be rewritten in the form of the Yule-Walker equations asap=−Rp−1rx,   (4)whereRp=XpXpH   (5)is the correlation matrix, andrx=XpxH   (6)is the correlation vector.
The relationship in Equation 4 can be formulated as an augmented Wiener-Hopf equation via a couple of algebraic manipulations. First, multiply both sides of Equation 4 on the left by Rp, then add rx to both sides to obtain
                                          [                                          r                x                            ⁢                              :                            ⁢                              R                p                                      ]                    ⁢                                          [                                                    1                                                                                      a                  p                                                              ]                =                              [            0            ]                    .                                    (        7        )            Adding the mean squared prediction error∥ε∥2=σx2+rxHap   (8)to the top row of Equation 7 yields
                                          [                                                                                σ                    x                    2                                                                                        r                    x                    H                                                                                                                    r                    x                                                                                        R                    p                                                                        ]                    ⁢                                          [                                                    1                                                                                      a                  p                                                              ]                =                              [                                                                                                                            ɛ                                                              2                                                                                                0                                                      ]                    .                                    (        9        )            Let ap+1 denote the augmented weight vector ap+1=[1 apT]T. Solving Equation 9 for ap+1 yields
                                          a                          p              +              1                                =                                    [                                                                    1                                                                                                              a                      p                                                                                  ]                        =                                                                              ɛ                                                  2                            ⁢                              R                                  p                  +                  1                                                  -                  1                                            ⁢                              u                                  p                  +                  1                                                                    ,                            (        10        )            where
                                          R                          p              +              1                                =                      [                                                                                σ                    x                    2                                                                                        r                    x                    H                                                                                                                    r                    x                                                                                        R                    p                                                                        ]                          ,                            (        11        )            and up+1, is a (p+1)×1 unit vector whose first element is unity with the remaining elements zero. The correlation matrix is now of dimension (p+1)×(p+1) and represents all the array elements. Once ap+1 has been found, it can be used in the power spectrum equation or power spectral density of an AR process:
                                                        P              x                        ⁡                          (              θ              )                                =                                                                    ɛ                                            2                                                                                                            ⅇ                    θ                    H                                    ⁢                                      a                                          p                      +                      1                                                                                                  2                                      ,                            (        12        )            where eθ is a steering vector directed toward the angle θ. For example, this steering vector can be a Fourier beamforming vector.eθ=einπ sin θ0≦n≦p   (13)
In most applications of the full rank AR just discussed, there is a compromise concerning the order of prediction, i.e., the number of samples to be incorporated into a linear combination and hence the number of unknowns in the normal equations or the number of parameters needed to characterize the all-pole model. A larger p provides a more accurate spectral envelope, with the disadvantage of requiring more computation, resulting in greater delays, and where the signal is transmitted, more bits for the transmission of the prediction coefficients. There is also a risk of over-modeling the process, which can lead to false peaks in the power spectral density. On the other hand, selecting a value for p that is too small results in no correct peaks.
Reduced rank AR estimation is known to alleviate a number of the disadvantages of full rank AR methods. First, reduced rank AR estimation can, in some cases, avoid over-modeling a process which is known to be of reduced rank. Second, the number of snapshots required to estimate the process is reduced. A process residing in a low dimensional space needs fewer measurements to form a good estimate. Third, reduced rank AR methods are less computationally demanding than full rank methods. The solution to AR spectrum estimation essentially involves inverting a pth order correlation matrix. If the data can be compressed into an M-dimensional subspace, inverting the resulting lower dimensional covariance matrix requires less computation.
To provide an example, if only two signals are present, a three element array should be able to detect both signals, albeit at a very low resolution. Each additional element of an array increases the resolution but may also give an additional false noise peak. Reduced rank techniques ideally select a subspace that captures signal information while excluding a large portion of the noise. Such methods require a lower level of sample support because they estimate the data in a smaller subspace. Algorithms which require fewer samples to estimate the process's statistics are particularly helpful when working in non-stationary environments.
Reduced rank processing may be thought as projecting the data matrix Xp into a reduced rank data space. The reduced rank data matrix isDM=LMLMHXp,   (14)where LM represents a matrix whose columns form an orthogonal basis that span the desired reduced subspace.
Using this new subspace data, the reduced rank AR weight vector can be calculated asap=−(DMDMH)†DMXH   (15)where † denotes the pseudo inverse.
Both the p-dimensional ãM vector and the M-dimensional ãM vector are termed “reduced rank vectors” since both were created in a subspace of rank M. Therefore, reduced rank AR spectrum estimation is the method of forming the best estimate of a pth order AR process in an M-dimensional subspace.
Two types of reduced rank AR estimators based on the principal eigenvectors of the data correlation matrix, i.e., principal component (PC)-AR, have been proposed in an effort to achieve the increased resolution of a high order full rank (FR)-AR model along with the fewer spurious peaks of the lower order FR-AR models. The first PC-AR method is a signal independent method which finds the eigenvectors of the Rp+1 correlation matrix and applies a reduced rank version of Equation 10. The second PC-AR method is signal dependent, which decomposes Rp into its eigenvectors and uses a reduced rank version of Equation 4. Both PC-AR methods take advantage of the diagonalized covariance matrix found in an eigenvector subspace so that no additional transformations are necessary to compute the inverse.
A fundamental limitation of the reduced rank PC-AR techniques is the need to work in a subspace having a rank that exceeds the unknown signal rank. If the subspace rank selected is less than the number of spatial frequencies present, the PC-AR method provides very poor performance. Thus, one must either have advance knowledge of an upper limit of the number of frequencies present or select a relatively high signal subspace to ensure that its rank is greater than the unknown number of frequencies, such that the small degree of rank reduction provides little benefit. Further, the need to perform eigenvector decomposition, even in reduced rank, requires a significant amount of computation to perform the full covariance matrix eigenvector decomposition, introducing delays which, in some situations, may be critical where rapid results are required.
Accordingly, the need remains for a method for reduced rank AR estimation which does not suffer the performance drawbacks of existing reduced rank methods and which requires less computational resources and time. The system and method disclosed herein are directed to such a need.